Universal Theory of Light Scattering of Randomly Oriented Particles: A Fluctuational-Electrodynamics Approach for Light Transport Modeling in Disordered Nanostructures

Disordered nanostructures are commonly encountered in many nanophotonic systems, from colloid dispersions for sensing to heterostructured photocatalysts. Randomness, however, imposes severe challenges for nanophotonics modeling, often constrained by the irregular geometry of the scatterers involved or the stochastic nature of the problem itself. In this Article, we resolve this conundrum by presenting a universal theory of averaged light scattering of randomly oriented objects. Specifically, we derive expansion-basis-independent formulas of the orientation-and-polarization-averaged absorption cross section, scattering cross section, and asymmetry parameter, for single or a collection of objects of arbitrary shape. These three parameters can be directly integrated into traditional unpolarized radiative energy transfer modeling, enabling a practical tool to predict multiple scattering and light transport in disordered nanostructured materials. Notably, the formulas of average light scattering can be derived under the principles of fluctuational electrodynamics, allowing the analogous mathematical treatment to the methods used in thermal radiation, nonequilibrium electromagnetic forces, and other associated phenomena. The proposed modeling framework is validated against optical measurements of polymer composite films with metal-oxide microcrystals. Our work may contribute to a better understanding of light–matter interactions in disordered systems, such as plasmonics for sensing and photothermal therapy, photocatalysts for water splitting and CO2 dissociation, photonic glasses for artificial structural colors, and diffuse reflectors for radiative cooling, to name just a few.

1 Theory of average light scattering from random oriented particles 1

.1 Formulation for a single object
In the context of Lippmann-Schwinger approach, 1 the presence of a object in an infinitely extended medium 0 induces a perturbation to the incident (source) electric fields, E i , that results in a scattered (induced) electric field, E s , in medium 0. For simplicity, we consider S2 medium 0 as vacuum, with electrical permittivity, ε 0 , and magnetic permeability, µ 0 , while the object is characterized by linear electromagnetic properties, i.e., dielectric constant, ε, and relative permeability µ, which can be, in general, nonlocal complex tensors. The total electric field, E t , that results from the superposition of E i and E s , satisfies the following where k 0 = ω √ ε 0 µ 0 is the wavevector in free space, I is the identity operator, H 0 = ∇ × ∇×, is the induced potential by the object, with V = 0 anywhere outside the object.
Solution of Eq. (S1) is given by: where G 0 is the free space dyadic green function: 2 G 0 (r, r ) = I + 1 which satisfies [H 0 − I] G 0 = δ(r − r ), and T is the scattering operator defined by: 1,3 T = 0 anywhere outside of the object n.
We begin our discussion with the definition of the absorption, P abs , and scattering, P sca , S3 energy fluxes in terms of the T-operator: 4 whose derivation stems, respectively, from the work done by the total (E t ) and scattered (E s ) fields over the polarization current, J s . Specifically, we consider the following relations: 4 The absorption and scattering cross sections [Main text, Eqs. (1a) and (1b), respectively] are given respectively by, 5 C abs = P abs / |E 0 | 2 2Z 0 , C sca = P sca / |E 0 | 2 2Z 0 . As discussed in the main text, to compute C abs and C sca , we need an expression for For an incident field of the form E i = E 0 e ik 0k i (r−r 0 ) , where r 0 represent the origin of the field away from the object, and considering the two orthogonal polarizationŝ andp: where the second step exploits the reciprocity of e ik 0k i ·(r−r ) , and the fourth step considers the identity, 6
The asymmetry parameter is given by the force from scattering fields, F s , through the relation: 7 The projection of the scattering force,k i · F s , expressed in terms of the induced currents is given by where we apply the relations from Eq. (S7).
Eq. (S10) represents a generalized form of the asymmetry parameter, regardless of the form of the source fields. Similarly to Eqs. (S6a) and (S6b), computation of the orientation average of,k i · F s , requires a direct expression for, k i where, in the third step, we consider the identity ∂ j G † 0 (r , r) = −∂ j G † 0 (r , r). 9 Substitution S5 of Eqs. (S11) and (S10) into Eq. (S9), together with the relation µ sca = µscaC abs Csca leads to the expression Eq. (4c) in the main text.
The expression Eqs. (4a), (4b) and (4c) in the Main text, can be extended for linear, homogeneous and non-absorbing external medium, by replacing k 0 , for n 0 k 0 , where n 0 is the refractive index of the external medium.

Surface-current BEM formulation
Particularly for the BEM, under the surface-current formulation, the fields and currents on an object n are, respectively, represented by the bi-linear expressions 10,11 where K and N are the surface electric and magnetic current, respectively.
Using bi-linear expressions, the electromagnetic field is given by φ n = Γ n ξ n . In this formula, where C n = i kn ∇ × G n , and k n and Z n are the wavevector and impedance of medium n, respectively.
Eq. (S22) in its bi-linear form is now expressed as: Notice that in this case, the elements of V and G 0 , are respectively given by V nm = δ nm Γ −1 n , and The solution is obtained by expanding the surface currents on a particular basis, f n r (r), as: ξ n (r) = r x n r f n r (r), where the expansion elements x n r are obtained from boundary conditions. Further details can be found elsewhere. 10,11 From the surface current expansion, the operators-based form of Eqs. (S25a), (S25b) and (S25c) is now given in matrix notations by: where the elements of the G 0 and G in are respectively given by G 0,nm In this notation, the operator , , denotes the conjugated inner product: where u and v are vector fields.

Branch No
(nm) (nm)  Figure S1: Average light scattering of agglomerated gold nanostars (a) C abs , C sca and µ sca of a single gold nanostar (grey) and cluster of agglomerated gold nanostars (red). The curves C abs and C sca are normalized to the number of stars (N star ) for direct comparison. The cluster is based on 15 stars at 20% volume fraction. (b) Model of the star considered for the simulations. The model consist on a star with 7 legs with the dimensions indicated by the table at the right-hand side of the figure. All the simulations were performed by the BEM code for average light scattering simulations. 12 The refractive index of gold can be found elsewhere. 13 Eqs. (S13a), (S13b) and (S13c) were implemented into a BEM application for average light scattering simulations. 12 The flexibility of the BEM algorithm enables to study objects of arbitrary morphology. 10 As a demonstration, we simulated the average light scattering parameters C abs , C sca and µ sca of a single gold nanostar and a cluster of agglomerated gold nanostars (Fig. S1). Fig. S1(a) illustrates the effects of interparticle coupling and interference of the scattered fields from the stars in cluster, which results in a reduction of C abs /N star and an enhancement of C sca /N star in comparison with a single star. Similarly, the scattering anisotropy is also clearly affected in the star cluster, as indicated by the variations in µ sca . Fig. S1(b) shows the dimensions of the gold nanostar model considered S8 for the simulations.

Average light scattering in Spherical wave basis: T-matrix formulation
The operator T for spherical wave basis is given by where P = M, N correspond to the two orthogonal polarizations, l = 1, 2...∞ and m = −l, ..., 0, ..., +l; T P P lm,l m are the elements of the T-matrix, 2 and E reg P lm (r) are the spherical waves regular at the origin, defined at the spherical coordinates r, θ and φ: 1 j l is the spherical Bessel function of order l, and Y m l (θ, φ) are the spherical harmonics according to the definition from Ref. [14]. Note that E reg * P lm (r) = E reg N l−m (r), where * is the complex conjugate. 1 Because G 0 is a symmetric and reciprocal dyadic, Asym (G 0 ) = Im (G 0 ). 9 For spherical wave basis, it can be demonstrated that: 1 Replacing Eqs. (S14) and (S15) into Eq. (4a) and (4b) (Main text), we derive the expres-S9 sions of C abs and C sca for T-matrix, i.e.: 15,16 C abs = 2π k 2 0 P lm P l m Re T P P lm,l m − |T P P lm,l m | 2 (S16) where in the C abs formula, we consider the relation To derive the T-Matrix formula of µ sca , we use the relation: which stems form the reciprocal and anti-symmetric properties of ∂ j G 0 . 9 For spherical waves basis, this expression is: 1 where p j P P lm, l m = ∂ j V P P lm,l m (a)| a=0 , and V is the translation operator of regular waves. 17 The form of V for spherical waves is given in Ref. [17], Appendix C.3.
By replacing Eqs. (S14) and (S18) into Eq. (4c) of the Main text, we derive the T-matrix formulation of µ sca : 0 j=x,y,z P lm P l m P l m p j P P lm, l m T * P P l m ,l m p j P P l m , l m T P P l m ,lm (S19) For a spherical object, Eq. (S19) can be simplified to: 1,18 (S20) S10 whereP = N if P = M and vice versa, T P l are the mie-scattering coefficients, 5 and: Using the identities

Formulation for multiple objects
In the case of multiple objects Eq. (S1), becomes: where N is the total number of objects. The solution of this equation is given by: 1 E t = E i + n G 0 T n E i n , and in matrix form: where φ i is a column vector, whose elements contain the incident fields on each object, is the analogous of T for a cluster; V is a band matrix operator whose elements are V nm = V m δ nm ; and the matrix operator G 0 represents the interaction between the objects in the free space, whose elements are G 0 nm = G 0 (ω; r n , r m ). The rest of the elements are defined in the main text.
Following the matrix notation, we extend the relations from Eq. (S7) to a vector form in S11 terms of the induced current, ξ s , and fields, φ s : Analogous to the derivation of C abs , C sca and µ sca for a single object, we use Eq. (S23) and the relations in Eq. (S24), to derive: these expressions represent the most general form of the average light scattering theory.
They can be applied to single and group of particles, including heterogeneous collections of particles.

Formulation for an individual object in a cluster
Using Eqs. (S25a), (S25b) and (S25c), we deduce the individual contribution of an object n in the cluster: where S = W Asym (G 0 ) W † , and S ∇ j = W Sym (∂ j G 0 ) W † .
In Eqs. (S26b) and (S26c) the first and second term inside the curly brackets represent, self and interference scattering, respectively. 19 Note that µ n sca in Eq. (S26c) is scaled by N S12 for better comparison with µ sca from an isolated object.

Approximation for subwavelength particles
For small objects, Asym (G 0 ) ≈ k 0 6π I, 20 and, T ≈ 4πk 0 α α α, a where α α α is the polarizability tensor. Substitution into Eq. (S26a), together with the relation gives: Because the particles are small, the second term in Eq. (S27) is negligible in comparison with the first, leading to C abs ≈ 1 3 (C abs,x + C abs,y + C abs,z ) .
The second term in Eq. (S27) corresponds to the scattering of the particle, which can be written as: where α ij are the elements of the tensor α α α. If α α α is a diagonal tensor, Tr α α α † α α α = |α xx | 2 + a the expression is deducted from the electric field generated by a point dipole 14 S13 |α yy | 2 + |α zz | 2 , and we derive the commonly used relation, 21,22 C sca ≈ 1 3 (C sca,x + C sca,y + C sca,z ) . and (4b) (Main text), respectively, and was computed by our BEM application for average light scattering simulations. 12 The approximated results are based on the relations, C abs ≈ 1 3 (C abs,x + C abs,y + C abs,z ) and C sca ≈ 1 3 (C sca,x + C sca,y + C sca,z ). The inset shows the relative error between the approximated and the exact solutions for small size parameters. The dimensions of the ellipsoid are indicated at the right figure. The refractive index of the ellipsoid is N = 3.5 + 0.1i.
In Fig. S2, we compute C abs and C sca for a randomly oriented ellipsoid, using the exact solution [Main text, Eqs. (4a) and (4b), respectively] and the approximation for small objects. The ellipsoid represents a typical case where α α α is not diagonal. As shown in the inset of Fig. S2, the approximation C sca ≈ 1 3 (C sca,x + C sca,y + C sca,z ), induces considerable error, even when k 0 L c is small. On the other hand, the approximation for C abs , shows a small relative error (< 2%) for small particles (k 0 L c < 2). For larger particles (k 0 L c > 2), both approximations fail. S14 dynamics 2.1 Vacuum friction and averaged light scattering An isolated object moving at a small velocity v experiences a non-conservative friction force, F f , that results from the interaction with thermal fluctuations in the free space. To a first order approximation, the vacuum friction is given by, F f = −γ · v, whereγ is the friction tensor given by: 18 here, k B is the Boltzmann constant, and T is the equilibrium temperature of the system.
Integrating the vacuum friction force over all solid angles leads to the relation, 23 where γ = iγ ii is the friction coefficient.
From the relations, ∂ i Asym (G 0 ) = −iSym (∂ i G 0 ), 18 and, the friction coefficient is expressed as: where F represent the number microcrystals of a particular morphology in the ensemble, and the subscripts "W, L" and "flake", indicate a microbars of dimensions L and W , and flakes, respectively. For F W,L , we considered the estimated size distribution [Main text, Fig. 3(c)], while F flake = 5, as an approximate of the number of flakes found in the sample.

S18
The volume of the microcrystal ensemble V is given by: where V W,L is the volume of a single bar with dimensions W and L, and V flake is the volume of the flake structure.

Prediction of refractive index and scattering from PE films
Overall, the optical properties of PE are highly dependent on the crystallinity and the manufacturing process. 27 Thus, we extracted the complex refractive index of a pure PE film (101 µm thick), n PE , fabricated under the same conditions than those used for the composite films. First [ Fig. S6(a)], T tot , T spec and R tot of the PE film where measured by FTIR and a gold integrating sphere (Main text, Materials and Methods). The real part of n PE was obtained from Palik et al. 27 The extinction coefficient, κ PE = Im (n PE ), was extracted from the equation: 28 where t film is the thickness of the PE film and R f is the reflectance at the air/PE interface.
R f was estimated from Fresnel law and the refractive index of PE from Palik et al. 27 As shown in Fig. S6(a), T tot = T spec at most parts of the spectrum, indicating a small scattering component in the film. To include this features in the modeling, we assume a small concentration of scattering particles inside the film. Because κ PE Re (n PE ), 27 Additionally, at low particle concentrations, 29 Thus, the scattering properties of the particles are extracted from: The extracted values of C sca for f v = 0.1% v/v and particles of 1 µm diameter are shown at the inset of Fig. S6(b). We fitted the results with a curve of the form C sca = a/λ b , with a = 2.4782 and b = 1.4095 [ Fig. S6(b)]. Finally, using Monte-Carlo simulations (see Main Text, Materials and Methods) we iterate to find the value of µ sca that best matches the results of experiments. We found excellent agreement for µ sca = 0.75 [ Fig. S6(a)].

Validation of Monte-Carlo Code
We validated our Monte-Carlo code for radiative transfer simulations against the Addingdoubling method [ Fig. S8(b)]. We calculated the total transmittance (T tot ) and reflectance (R tot ), and the specular transmittance (T spec ) at normal incidence over a 1 mm thick film with refractive index n film = 1.4, and 0.1% particles per volume. The particles have a diameter D p = 1 µm, and the scattering properties shown in Fig. S8(a). S21 a b Figure S8: Validation of Monte-Carlo code for radiative transfer calculations against Adding-doubling method. (a) Light scattering properties C abs , C sca and µ sca of the spherical particles (1 µm diameter) considered in this study (b) Total transmittance and reflectance and specular transmittance (inset) for light at normal incidence of a 1 mm thick film with 1.4 refractive index and 0.1% particles per volume. The results from our Monte-Carlo code (circles) are plotted against Adding-doubling method calculations using the code by Prahl. 30 The results for Adding-doubling method are obtained from the open-source code by Prahl. 30 This code was used to calculate the total reflectance and transmittance. The specular transmittance is calculated as where R = 4.12% is the reflectance of a 1 mm thick film with 1.4 refractive index, f v is the volume fraction, V p is the particle's volume, and t film is the thickness of the film.
The results of total transmittance/reflectance from our Monte-Carlo code show excellent agreement with the curves from Adding-doubling method. The specular transmittance (inset of Fig. S8) shows good agreement up to 10 nm wavelength. Outside this range, the values of specular transmittance are bellow the minimum resolution of our Monte-Carlo setup (0.0001% for 1,000,000 photons).